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ABSTRACT 

We explore the dynamics of a "cluster wind" flow in the regime in which the shocks resulting from 
the interaction of winds from nearby stars are radiative. We first show that for a cluster with T Tauri 
stars and/or Herbig Ae/Be stars, the wind interactions are indeed likely to be radiative. We then 
compute a set of four, three dimensional, radiative simulations of a cluster of 75 young stars, exploring 
the effects of varying the wind parameters and the density of the initial ISM that permeates the volume 
of the cluster. These simulations show that the ISM is compressed by the action of the winds into 
a structure of dense knots and filaments. The structures that are produced resemble in a qualitative 
way the observations of the IRAS 18511+0146 of Vig et al. (2007). 
Subject headings: Hydrodynamics - shock waves - stars: winds, outflows 



1. INTRODUCTION 

In a recent paper, Vig et al. (2007) reported JCMT- 
SCUBA and Spitzer (IRAC & MIPS) observations of 
the region around IRAS 18511+0146, and deduced that 
there is a young cluster of Herbig Ae/Be stars (probably 
also having lower mass stars), which is still embedded 
in a dense, highly extinguished region. The interstellar 
medium (ISM) within this cluster has a complex struc- 
ture of dense clumps and filaments, which are seen as 
structured extinction of the background infrared (IR) 
emission. In the present paper, we explore whether or 
not this clump/filament structure could be the result of 
the stirring up of the ISM by the winds from the young 
stars. 

This cluster appears to be a younger, more embedded 
version of clusters of Herbig Ae/Be stars such as the ones 
described by Testi et al. (1997, 1999). These clusters 
typically have less than 100 stars within regions of ~ 1 pc 
size (corresponding to I) ~ 0.1 pc separations between 
the cluster stars). 

The region between the stars in this kind of clus- 
ter will of course be stirred up by the massive winds 
ejected by the young stars. These winds will have high, 
> 10~^ M0yr~^ mass loss rates and terminal veloci- 
ties = 100 — 400 km s~^. These parameters, together 
with the D ^ 0.1 pc separations between the cluster stars 
(see above) imply that the shock interactions between the 
winds from nearby stars might be radiative. 

The "cluster wind" resulting from the combined effect 
of all of the stellar winds will in this case be very different 
from the non-radiative cluster winds studied by Chevalier 
& Clegg (1985) and Canto et al. (2000). These papers 
have all studied cluster winds with no radiative losses 
(Chevalier & Clegg 1985, Canto et al. 2000, Raga et al. 
2001 and Rockefeller et al. 2005, Rodriguez- Gonzalez et 
al. 2007), or with relatively small radiative losses (Silich 
et al. 2004, Tenorio-Tagle et al. 2005), and are applica- 
ble for clusters of massive stars (with wind velocities of 
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~ 1000 km s ^ which result in non-radiative wind- wind 
interactions). 

Motivated by the observations of IRAS 18511+0146 
of Vig et al. (2007), we have studied the flow result- 
ing from a cluster of stars with wind interactions that 
produce radiative shocks. As this problem has not been 
studied before, we present a series of idealized models 
in which all of the stars in the cluster have identical, 
isotropic winds, and in which the initial ISM permeat- 
ing the volume of the cluster is homogeneous. Because 
of these simplifications, our models should be regarded 
as an exploration of the properties of highly radiative 
cluster wind flows rather than an attempt to model a 
particular object (namely, IRAS 18511+0146). 

In the present paper, we first study for which combina- 
tions of parameters (mass loss rate M^, stellar wind ve- 
locity Vw , and separation D between nearby cluster stars) 
one obtains a highly radiative cluster wind flow (§2). We 
then compute four 3D simulations of cluster winds in 
the highly radiative regime (§3), and obtain predictions 
of the highly structured column density maps that result 
from our simulation (§4). Finally, in §5 we present our 
conclusions. 

2. RADIATIVE LOSSES IN A CLUSTER WIND 

Let us consider a stellar cluster with a local stellar 
density n (=number of stars per unit volume), of stars 
with identical, isotropic winds with a mass loss rate M^j 
and a terminal wind velocity Vw The typical separation 
between stars then is D = n~^/^. 

The two-wind shock interactions between nearby stars 
occur at a typical distance ^ D/2 from each of the stars, 
so that the typical pre-shock densities have values 

_ 

where uih is the H mass, and we have assumed a 90% H 
and 10% He particle abundance. 

The shock interactions between nearby stars will be ra- 
diative if the cooling distance dcooi satisfies the condition 
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In order to estimate the cooling distance, we proceed 
as follows. We consider the head of the stellar bow shock 
structures that will be formed with a shock velocity of 
^ Vyj (i. e., equal to the stellar wind velocity) and a 
preshock density given by equation Q. We then use the 
fact that the cooling distance (for shocks with cooling 
in the low density regime) scales as the inverse of the 
pre-shock density to write 



dcool ij^pre 1 ^) 



100 cm-^ 



^c,100 



iv): 
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where v {= v^) is the shock velocity, and dc^ioo{v) is the 
cooling distance behind a shock with Upre = 100 cm~^. 
We then estimate the function (ic,ioo('^) by carrying out 
a fit to the cooling distance (to 10^ K) obtained from 
the "self-consistent preionization" , steady, plane-parallel 
models of Hartigan, Raymond & Hartmann (1987). We 



use 
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which fits the computed cooling distances with a rela- 
tive accuracy of better than 15% in the v = 100 
400 km s~^ shock velocity range. We note that large 
deviations between this fit and the values of Hartigan 
et al. (1987) occur for v < 100 km s~^, so that the fit 
should not be applied for lower shock velocities. 



Combining equations 
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In Figure 1, we plot the cooling parameter n (see equa- 
tion [s]) as a function of Vy^ and My^/D. This figure shows 
that a substantial part of the parameter range that would 
be expected for a cluster of low or intermediate mass stars 
produces cooling parameters n <1. 

A cluster wind in this "high cooling regime" will have 
dense, cool structures resulting from the radiative shocks 
in the interactions between nearby stars. In the following 
section, we present a numerical simulation of a flow in 
this regime. 

3. A 3D SIMULATION OF A RADIATIVE CLUSTER 

WIND 

3.1. The numerical setup 

We have carried out 4 numerical simulations in which 
we place 75 stars which are uniformly distributed inside 
a sphere of radius Rc = 0.5 pc. The stars have identical 
mass loss rates = 10"^ (models Ml and M2), 4x10"^ 
(M3) or 10"^ Moyr-^ (M4) and wind velocities Vu, = 100 
(Ml), 200 (M2 and M3) or 300 km s'^ (M4). These 
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Fig. 1. — This diagram shows the "coohng parameter" hz = 
dcool/D (where D is the typical the separation between stars) as a 
function of Mw x (O.lpc/D) and Vw, where Mw is the mass loss rate 
and Vw the wind velocity. The thick curve shows the = 1 contour, 
the other contours represent values of k at successive factors of two 
below (top left region of the diagram) or above (bottom right region 
of the diagram) unity. The region above the thick curve therefore 
represents the parameter space in which the two-wind interactions 
in a cluster wind flow are radiative. 



parameters are listed in Table 1. These wind parameters 
span a parameter range which is appropriate for low and 
intermediate mass stars, and have been chosen so as to 
produce a cluster wind in the highly radiative regime (see 
below). 

With the number of stars and the volume of the clus- 
ter, we obtain a density of stars n = 133.7 pc~^, thus a 
typical separation between stars of D = n~^^^ = 0.20_pc. 
For these values for and v^j^ using equations ([5|6| 

we obtain k, = 0.008, 0.518, 0.130 and 0.666 for moaels 
M1-M4, respectively. In other words, the post-wind bow 
shock cooling distances have values ~ 0.5 ^ 30 % of the 
typical separations between stars, placing the simulations 
in the highly radiative regime. 

The stars are placed randomly within the volume of the 
spherical clusters, elliminating the positions that appear 
within a distance of < 0.035 pc from any of the previ- 
ously generated stellar positions. We then impose the 
isotropic stellar wind condition within spheres of radii 
Ryj = 0.0172 pc centered at each of the stellar posi- 
tions. A uniform, = 5000 K temperature is imposed 
within each of the wind sources. The rest of the com- 
putational domain is filled with a homogeneous medium 
with number density Uenv = 10^ (models Ml and M2) or 
5x 10"^ cm" ' 
perature Te, 

simultaneously at the beginning of the time-integration. 
The winds and the environment have neutral H, and a 
seed electron density which is assumed to come from 
singly ionised Carbon. 

The gasdynamic equations, together with a rate equa- 
tion for neutral Hydrogen are integrated using the 



^env 

(models M3 and M4, see Table 1) and a tem- 
Lenv = 100 K. The stellar winds are "turned on" 
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"yguazii-a" code (Raga et al. 2000, 2002). The energy 
equation includes the cooling function described by Raga 
& Reipurth (2004). This cooling function is appropriate 
for describing the cooling of the shocked wind material. 
However, it is not necessarily appropriate for the shocked 
ISM, which is likely to be initially molecular. At the rela- 
tively low resolution of our simulations, however, the de- 
tails of the cooling function are unlikely to be important, 
as the cooling regions are at best marginally resolved. 

The Tenv = 100 K chosen for the initial configuration of 
the environment is higher than the ~ 10 K temperature 
expected for a molecular cloud. This artificially high 
choice of temperature is consistent with the fact that 
our cooling function (see above) is appropriate only for 
partially ionized gas, as it does not include the cooling 
processes important at temperatures ~ 10 K. 

The numerical integrations are done in a cubic domain 
with a size of 1.2 pc, using a 5-level, binary adaptive grid 
with a maximum resolution of 512 (Ml) or 256 (M2-M4) 
points along each of the three axes. The maximum res- 
olution level is only allowed within the spheres in which 
the stellar wind conditions are imposed. Outfiow condi- 
tions are imposed on all of the grid boundaries. 

Together with the gasdynamic equations and the rate 
equation for neutral H (see above), we have integrated 
an advection equation for a passive scalar. This scalar 
has a positive value for the environmental material, and 
a negative value for the stellar winds. In this way, we can 
at all times trace which regions of the computational do- 
main are occupied by wind or by environmental material 
(identified by the sign of the passive scalar). 

3.2. Model results 

We have carried out time integrations from t = (when 
the stellar winds are "turned on") up to t = 1.9 x 10^ yr 
for models M1-M4. Figure 2 shows a "volume-rendered" 
depiction of the 3D density structure of model Ml at this 
final time. We see that there is a dense shell (which has 
partially left the computational domain) which is being 
pushed out by the stellar winds into the surrounding en- 
vironment. This shell is formed by material from the 
winds and by part of the ISM material which was ini- 
tially permeating the volume of the cluster. Also, we see 
that there is a complex structure of filaments and clumps, 
which are formed from material that was initially present 
in intercluster medium and wind material that has gone 
through the stellar wind bow shocks and cooled to low 
temperatures (and high densities). 



TABLE 1 
Model parameters 



Model 






^cloud 




resolution 




km s~-^ 


M© yr-i 


crQ— 




pixels 


Ml 


100 


10-6 


10^ 


0.004 


512^ 


M2 


200 


10-6 


10^ 


0.259 


256^ 


M3 


200 


4x10-6 


5x10^ 


0.065 


256^ 


M4 


300 


10-5 


5x10^ 


0.333 


256^ 



^Using the fit to the models of Hartigan et al. (1987), i.e. 
equation (Is} 



In Figure 3, we show column density maps (obtained 
by integrating the number density along the z-axis of the 
computational domain) at integration times t = 19, 3.2, 



11 and 3.2 x 10"^ yr for models Ml, M2, M3 and M4, 
respectively. These times have been chosen so that the 
outer, dense shell is starting to leave the computational 
domain in each of the models. All of the maps have a 
central condensation, which has a peak column density of 
9.11, 12.8, 66.6 and 36.3 x 10^^ cm-^ for models M1-M4, 
respectively. 

Figure 4 shows column density maps for models Ml- 
M4 computed only with the material from the ISM which 
originally permeated the computational domain. In these 
maps, the central condensation has a peak column den- 
sity of 1.82, 10.8, 29.7 and 33.6 x 10^^ cm-^ for models 
M1-M4, respectively. Comparing the peak column den- 
sities for ISM material only (Figure 4) and the column 
densities for ISM+wind material (Figure 3), we see that 
while for model Ml the main contribution to the col- 
umn density of the central condensation comes from the 
wind material, for models M2-M4 the main contribution 
to this column density comes from the material in the 
initial ISM. 

We have also computed the fraction fm of mass of the 
initial ISM which is still present within the volume of 
the cluster (i. e., within a sphere of radius Rc = 0.5 pc) 
for the time frames shown in Figure 4. We obtain 
fm = 0.076, 0.026, 0.26 and 0.052 for models M1-M4, 
respectively. Therefore, for model M3, almost 75% of 
the initial ISM mass has already been expelled from the 
volume of the cluster (most of this mass being present 
in the expanding shell pushed out by the cluster wind 
fiow), and for the other three models more than 90 % of 
the ISM mass has been expelled. 

As the material from the winds is unlikely to have dust, 
the column density distributions of the ISM material (see 
Figure 4) are proportional to the expected extinction. If 
we use a standard conversion factor Av = IO'^^Nh, 
where Ay is the optical extinction and Nh is the hy- 
drogen column density in cm~^, the peak values of Ay 
corresponding to the central condensations are Ay ^ 18, 
110, 300 and 340, for models M1-M4, respectively. 

Finally, in Figure 5 we present a series of column den- 
sity maps (computed only for the material of the initial 
ISM) illustrating the time-evolution obtained from model 
M3. From this time-sequence, we see that we can already 
identify the structure that forms the central condensa- 
tion at t = 1.6 X 10^ yr. This condensation persists for 
~ 10^ yr (i. e., to the end of the time-integration). Sev- 
eral other features can be seen in both the t = 4.8 x 10^ 
and t = 7.9 x 10^ yr frames, indicating that the filaments 
have lifetimes ~ 10^ yr. 

4. CONCLUSIONS 

Recent Spitzer observations of a young cluster of Her- 
big Ae/Be stars (Vig et al. 2007) show that the IR extinc- 
tion towards IRAS 18511+0146 has a curious structure 
of filaments and clumps. This observation has motivated 
us to study the formation of dense, neutral structures 
as a result of the interaction between the winds of the 
young stars. 

We find that if the Herbig Ae/Be stars have wind veloc- 
ities < 200 km s~^, the stellar wind interactions be- 
tween nearby stars are highly radiative for mass loss rates 
as low as ^ 3x lO~^M0yr~^ (see Figure 1 and equa- 
tion[5|. However, for higher wind velocities, considerably 
higher mass loss rates are required for the wind interac- 
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The cubic domain has a size of 1.2 pc. 



tions to be radiative. For example, for = 400 km s ^, 
mass loss rates > lO~^M0yr~^ are required (see 
Figure 1). The values quoted in this paragraph are for 
an average separation of 0.1 pc between the stars in the 
cluster. 

These values indicate that in a cluster of Herbig Ae/Be 
stars many of the wind interactions between nearby stars 
are likely to be in the highly radiative regime. We then 
compute a 3D simulation of a cluster wind flow in this 
regime, and we show that it does lead to the production 
of a structure of dense filaments and clumps. 

From our simulations, we obtain predicted column den- 
sity maps which show spatial structures which are quali- 
tatively similar to the observations of IRAS 18511+0146 
(Vig et al. 2007). Both the predicted and the observed 
maps show a structure of dense clumps and filaments, 
distributed within as well as outside the volume of the 
clump. We find that in all simulations a central clump 
is produced, which has a lifetime of ~ 10^ yr. 

It is more difficult to carry out a more quantitative 



comparison between our predictions and the observations 
of IRAS 18511+0146. Vig et al. (2007) first suggest 
that the filament /clump structures in this object have 
Ay ~ 6-8 extinctions (deduced from their Spitzer IR 
observations). However, when they compare these val- 
ues with the extinction deduced from millimetre obser- 
vations, they conclude that the IR maps probably have a 
strong contribution from foreground emission, and that 
the real values of the extinction should be Ay ^ 50. 
Therefore, only our model Ml produces visual extinc- 
tions which are too low compared with the ones deduced 
from the observations. 

We end by noting that the simulation presented in this 
paper is the first attempt that has been made for mod- 
elling the interaction of the outflows from a cluster of 
young, low and intermediate mass stars. Our simulation 
was done assuming that : 

• the stellar winds are "turned on" simultaneously, 

• the region is initially filled by a homogeneous 
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Fig. 3. — Column density maps obtained from models M1-M4 by integrating the number density along the 2;-axis. The stratifications 
are shown with the logarithmic gray-scale (color-scale in the online version) given by the bar on the top right plot. The integration time 
for each model is indicated in the label at the top of each panel. 



medium, 

• the stellar outflows are isotropic, 

• all of the stellar winds are identical. 

It is particularly important to remove this last approxi- 
mation if one wants to model IRAS 18511+0146 in de- 
tail, as this cluster harbours the massive protostar IRAS 
18511 A. 

The other assumptions in our models (listed above) 
are also not realistic, and the effects of removing them 
should be explored in future work. For example, assum- 
ing that the winds are isotropic is somewhat dubious in 
the context of a cluster of young stars. The more massive 
stars in such a cluster will evolve faster from an early, 
collimated outflow phase into a later phase in which a 
more isotropic wind is produced. The lower mass stars 
in the cluster will evolve less rapidly, and will remain in 
a collimated outflow phase for a longer period. For ex- 
ample, for a cluster with an age of ~ 10^ yr one might 
expect to have a few more massive (~ lOM©) stars with 
isotropic winds, and a larger number of lower mass stars 



still producing jets. This combination of isotropic and 
collimated winds will result in flow configurations that 
might be applicable when more detailed observations of 
the wind interaction in clusters of young stars become 
available. 

Another point that would be worthwhile to explore is 
the effect of the cooling function. As we have described 
in §2, our simulations are computed with a parametrized 
cooling function appropriate for partially (or fully) ion- 
ized gas. If one included a description of the non- 
equilibrium chemistry, and used it to compute a molec- 
ular cooling function, one would obtain cooling to lower 
temperatures in the dense filaments of compressed cloud 
material. This would lead to the production of narrower 
and denser filaments. A calculation including these ef- 
fects should also have a substantially higher resolution in 
order to be able to resolve the associated cooling regions. 

The conclusion that we obtain from the work described 
above is that the filamentary structure observed in IRAS 
18511+0146 by Vig et al. (2007) might be the result of 
the surrounding ISM being compressed into dense sheets 
by the wind interactions between the cluster stars. How- 
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Fig. 4. — ISM column density maps obtained from models M1-M4 by integrating the number density of the material in the initial ISM 
along the z-axis. The stratifications are shown with the logarithmic gray-scale (color-scale in the online version) given by the bar on the 
top right plot. The integration time for each model is indicated in the label at the top of each panel. 



ever, it is of course clear that at least part of the observed 
structures could have been present initially in the dense 
cloud from which the cluster stars were formed. 
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Fig. 5. — Time-sequence showing the evolution of the ISM column density structure of model M3 (obtained by integrating the density 
of the material in the initial ISM along the z-axis). The stratifications are shown with the logarithmic gray-scale (color-scale in the online 
version) given by the bar on the top right plot. 



